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Approximate  Exit  Probabilities  for  a  Brownian 
Bridge  on  a  Short  Time  Interval,  and  Applications 


by 


H.  R.  Lerche  D.  Siegmund 

University  of  Freiburg  Stanford  University 

To  Henry  Daniels  on  his  75th  Birthday 


Summary. 

\  k 

v  Let  T  be  the  first  exit  time  of  Brownian  motion  W(f)  from  ayegion  72.  in  d-dimensional 

Euclidean  space  having  a  smooth  boundary.  Given  points  add  £1  in  7(2,  ordinary  and 

>< ; 

large  deviation  approximations  are  given  for  Pr{T  <  ff|W(0)  =  $>,  W(t)  =  &}  as  t  — *  0. 
Applications  are  given  to  hearing  the  shape  of  a  drum,  approximating  the  second  virial 
coefficient,  and  Monte  Carlo  estimation  of  first  passage  distributions  fcir  Brownian  motion. 


Kev  Words  and  Phrases:  Brownian  bridge,  first  passage,  hearing  the  shape  of  a  drum, 
Monte  Carlo  methods. 


1.  Introduction. 


Let  W(f),  0  <  t  <  oo,  denote  Brownian  motion  in  lRd  with  W(0)  =  fo-  For  t  >  0  and 
events  A  in  the  <r-algebra  generated  by  W(.s),0  <  s  <  t,  let 


=  p^\W(o)  =  b,w<t)  =  6). 


Assume  that  £o  and  belong  to  some  region  72  with  a  smooth  boundary  #72,  and  let  T 
denote  the  time  W  first  leaves  72,  i.e.,  T  =  inf{t :  W(f  )e<572}.  The  principal  subject  of  this 
paper  is  the  asymptotic  behavior  of 


4Vr<«> 


as  t  — *  0  and  the  &  are  at  a  distance  0(t1/2)  from  each  other  and  from  #72.  A  secondary 
consideration  is  the  case  where  the  distances  of  the  &  to  the  boundary  and  each  other  are 


fixed  as  t 


This  problem  for  d  =  2  and  £o  =  £i  arises  naturally  in  the  beautiful  paper  of  Kac 
(1966),  who  was  concerned  with  the  behavior  for  small  t  of 


(1.2) 


Eexp(-Afet), 


where  the  A*  are  eigenvalues  of  the  Laplacian  acting  on  functions  having  domain  72  and 
vanishing  on  <572.  Kac  shows  that  as  t  — ►  0  (1.2)  has  an  expansion  of  the  form  c\t~l  + 
cit~xl2  +  o(<-1/2),  where  c\  and  c-i  are  numerical  multiples  of  |72|,  the  area  of  72,  and  |#72|, 
the  length  of  <572,  respectively.  He  argues  heuristically  that  the  next  term  is  (1— h)/6,  where 
h  is  the  number  of  holes  in  72.  For  more  detailed  results  along  these  lines,  see  Louchard 
(1968),  McKean  and  Singer  (1967),  Stewartson  and  Waechter  (1971),  and  Smith  (1981). 
Of  these,  only  Louchard  attempts  a  probabilistic  analysis,  and  his  argument  appears  to 


contain  a  mistake. 


Starting  from  the  physical  problem  of  evaluating  the  second  virial  coefficient  of  a 
hard  sphere  gas,  Handelsman  and  Keller  (1966)  arrive  at  essentially  Kac’s  mathematical 
problem,  for  the  case  d  =  3,£o  =  £i ,  and  72  the  region  exterior  to  a  sphere.  They  derive 
what  in  Kac’s  problem  corresponds  to  C3  and  the  next  term,  erf1/2 .  Although  their  method 


.  j  'y  'j'-  /  v.  ■/.  f.  V.  «\  -r.-r.  <r 


does  not  seem  capable  of  being  turned  into  a  rigorous  proof,  minor  modifications  appear 
to  produce  correct  answers  under  much  more  general  conditions. 

A  problem  having  a  rather  different  flavor  is  to  estimate  (1.1)  by  a  Monte  Carlo  exper¬ 
iment,  when  t  is  not  necessarily  small  and  the  £,  are  not  necessarily  close  to  the  boundary 
of  H.  A  natural  approach  is  to  partition  the  time  interval  [0,  f]  at  m  +  1  equally  spaced 
points  t,  —  it/m ,  t  =  0, . . . ,  m  and  count  the  relative  frequency  with  which  a  simulated 
path  W(ti),i  =  0,1,..., m,  leaves  H.  The  bias  introduced  by  discretization  is  typically 
0(m-1/2)  (cf.  Siegmund,  1985,  Chapter  X,  or  Hogan,  1984);  and  increasing  m  sufficiently 
to  reduce  this  bias  to  an  acceptable  level  is  computationally  time  consuming.  However, 
having  observed  TV (#,_! )  =  and  W(t,)  =  £,e7l,  one  can  use  an  approximation  to 

(1.1)  and  a  single  uniform  random  variable  to  simulate  the  event  that  W(s)  leaves  1Z  for 
some  s  in  the  time  interval  (<<_!,<,).  Although  the  original  interval  [0,t]  need  not  be  short, 
the  various  subintervals  are,  provided  m  is  large. 

Note  that  this  technique  does  not  require  that  W  be  exactly  Brownian,  but  only  that 
it  be  approximately  so  over  short  time  intervals.  The  basic  idea  is  in  principle  applicable 
to  diffusion  processes  and  to  certain  Gaussian  processes  which  are  locally  Brownian. 

As  noted  by  Kac,  in  the  case  £o  =  £i,  the  probability  (1.1)  is  to  a  first  order  approx¬ 
imation  equal  to  the  probability  that  W(a)  for  some  0  <  s  <  t  touches  the  plane  tangent 
to  V.  at  the  point  of  &R.  closest  to  £0-  Section  2  contains  the  first  term  of  an  Edgeworth 
type  expansion  for  this  probability  when  £o  and  are  not  necessarily  the  same.  A  large 
deviation  approximation  is  also  given.  Section  3  gives  the  substantially  more  complicated 
second  Edgeworth  term.  For  computational  simplicity  only  the  case  £o  =  £i  is  considered 
there,  but  this  case  illustrates  the  method  and  contributes  to  the  Kac  and  Handelsman- 
Keller  problems.  The  method  used  in  Sections  2  and  3  is  a  modification  of  that  introduced 
by  Siegmund  and  Yuh  (1982)  in  a  simple  linear  case  and  explored  more  thoroughly  by 
Siegmund  ( 1985 ) . 

Section  4  describes  some  illustrative  Monte  Carlo  experiments.  It  can  be  read  in¬ 
dependently  of  Sections  2  and  3,  except  for  an  occasional  reference  to  some  of  the  basic 
notation  and  to  the  statements  of  Theorems  1  and  2. 
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2.  Approximations  to  P^Jtil  {T  <  <}. 

For  ease  of  exposition  we  consider  in  detail  only  the  case  d  =  2  and  indicate  in  some 
remarks  how  the  results  are  modified  for  general  d. 

Given  £o>  £i<72,  close  to  &R,  and  to  each  other,  assume  there  exists  a  unique  point  on 
dTt  the  sum  of  whose  distances  to  the  given  points  is  a  minimum.  Consider  the  Cartesian 
coordinate  system  which  has  this  point  as  its  origin,  the  x-  axis  as  tangent  and  the  y-axis 
as  outward  normal  to  72.  There  exists  a  function  y  =  /(x)  such  that  locally  near  (0, 0)  <972 
is  given  by  the  graph  ( x,/(x )).  Let  £,  have  coordinates  (x*,  y,)  in  this  coordinate  system, 
and  assume  that  y*  <  0(*  =  0, 1).  It  is  easily  seen  that  and  satisfy  — xo/|yo|  =  xi  / |yi  |. 
(A  ray  of  light  emanating  from  £o  and  reflecting  off  the  x-axis  at  the  origin  passes  through 
.)  Let  W{t)  denote  Brownian  motion  starting  from  W(0)  =  £0>  and  define 

(2.1)  T  =  inf{t  :  W2(t)  >  /(Wi  (<))}. 

In  general,  T  is  not  the  exit  time  of  W  from  72,  but  for  £o  close  to  (0, 0)  it  is  with  probability 
close  to  one.  (A  more  precise  estimate  is  given  below.) 

In  order  to  study  ^  {T  <  c}  it  is  convenient  to  use  Brownian  scaling  to  replace  the 
given  problem  on  the  time  interval  [0,  e]  by  an  equivalent  one  on  [0,1].  Since  W(et)/e1/2 
is  Brownian  motion  starting  from  £o  =  £o it  is  easy  to  see  that 

(2.2)  />,';»  i(T<t)=P»'(-(f<l), 
where  =  ii/e1!2  (j  =  o,  1)  and 

(2.3)  T  =  f,  =  inf{<  :  W2(t)  >  t~xl2 

To  give  a  precise  statement  of  our  first  result  it  is  convenient  to  change  our  viewpoint 
slightly  and  regard  /  as  given  and  the  points  as  variable. 

Theorem  1.  Assume  /  is  twice  continuously  differentiable,  /( 0)  =  /'( 0)  =  0.  and 
/"( 0)  ^  0.  Suppose  =  (x,,y.)  (i  =  0, 1)  satisfy  y*  <  0  (t  =  0, 1)  and  —  x0/|yo I  =  J*i /|yi  U 
and  converge  to  (0,0)  as  e  — *  0  in  such  a  way  that  i,  =  i,/e x^2  are  fixed  (i  =  0. 1).  Then 
for  T  defined  by  (2.1) 
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(2.4)  P<',<,(T  <  s)  =«Jtp(-2yoyi/t)|l  -/"(°)  e  1/2|ws/i| 

%to+y,‘)/^!(i  -  £~i(i‘  - *•>’> + + °<£i/!)  }■ 
where  $  and  <p  Me  the  standard  normal  distribution  and  density  function  respectively. 

Proof.  By  virtue  of  (2.2)  it  suffices  to  consider  the  standardized  problem  on  the  time 
interval  [0,1],  with  fixed  initial  and  terminal  points  fo  =  (xo,yo)  and  =  (x^yi),  and 
T  defined  by  (2.3).  To  simplify  the  notation  we  consider  only  this  standardized  problem 
and  omit  the  tildes  for  the  rest  of  the  proof.  In  this  new  notation,  where  all  variables  have 
tildes,  but  the  tildes  are  omitted,  (2.4)  becomes 

(2.5)  {Tt  <  1}  =  exp(-2y0yi  )|  1  -  el/2f"(0)  |y0yi  I 

x  *;■  —  ~(xi  -x0)2)  +  (xg|yi|  +  xj|y0|)  +  o(l)  }, 

V?(yo  +  yi )  J  J 

with  Tt  defined  by  (2.3). 

We  begin  with  an  informal  calculation  and  provide  a  justification  later.  The  argument 
proceeds  from  a  suitable  likelihood  ratio  identity.  Let  =  (xj,  |yi|).  The  likelihood  ratio 
of  W(s), s  <  t,  imder  relative  to  P^\>  is  easily  calculated  to  be 

exp(-2|y0yi|)exp[-2jy1|W2(t)/(l  -  <)]. 

Thus  since  W2(T)  =  s~1^2  /(e1^2  Wj(T)),  we  have 


<  1}exP(2lyoyil) 


-  {' 


— 2|yi  \f{elf2W\{T)) 
e1/2(l  -  T) 


T  <  1 


(cf.  Siegmund,  1985,  Proposition  3.12). 
Since 


(2.7) 


e~'/2f{elf2x)  ~  c1/2/"( 0)x2/2  0  (e  —  0), 
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for  all  sufficiently  small  e,  {T  <  1}  =  1,  and  the  right  hand  side  of  (2.6)  can  be 
expanded  to  become 

(2.8)  1  -  In  |£'/V"(0)E<;;{,  [W?(T)/(1  -  T)]  +  . . .  . 

Define 


(2.9) 


r  =  inf{f :  Wi(t)  >  0}. 


From  (2.7)  follows  {Te  — ►  r}  =  1,  and  hence  (one  expects  that) 

(2.io)  \wut)k  i  -  t)]  -  WM/a  -  r». 

It  is  easy  to  see  that  conditional  on  W2(t),t  <  is  distributed  as  [r(l  —  r)]1/2Z  + 

x0  +  (x\  —  Xq  )r,  where  Z  has  a  standard  normal  distribution.  Hence 


(2.11) 


-  r>l  =  *o  +  U  +  2*o*i  -  *5)<VT) 

+  x?£<;;(;[rV(l-r)]. 


Equation  (2.5)  follows  from  (2.6),  (2.8),  (2.10),  (2.11),  and  the  evaluations  given  below  in 
Lemma  1. 

To  make  the  proceeding  manipulations  into  a  proof,  one  must  consider  the  remainders 
in  (2.7)  and  (2.8),  and  justify  the  convergence  indicated  in  (2.10). 

Let  A  =  {max«<r  |Wj(<)|  <  r-1/4}.  From  the  distribution  of  the  maximum  of  a 
pinned  Brownian  motion  (e.g.,  Siegmund,  1985,  (3.13)),  it  is  easy  to  see  that 


(2.12) 


=  °<£'>  for  *“  *  >  °- 


>0) 


Hence  (2.6)  can  be  replaced  by 


(2-13)  {T  <  1}  exp(— 2|y0yi|) 


-  *£«  { 


exp 


-2\yx\f{e'l*W,{T)) 
£1/2(1  -  T) 


;  {T  <  1}  n  A  >  +  o(sk) 


} 


for  all  k  >  0.  Let  6  >  0.  By  (2.12)  and  two  applications  of  Taylor’s  theorem  with  remainder 
along  the  lines  suggested  in  (2.7)  and  (2.8)  one  can  obtain  upper  and  lower  bounds  for  the 
right  hand  side  of  (2.13)  in  the  form 

1  -  Ivi  k1/2[/"(  0)  ±  «]£{”.  [Wl(T)/(\  —  T);  A]  +  o(ek). 

Since  8  >  0  is  arbitrary,  by  (2.11)  and  Lemma  1  below  it  suffices  to  show  (cf.  (2.10)) 

Wcn/d  -  n  4  -  Wdt/d  -  »•)], 

where  r  is  defined  by  (2.9).  Since  P^(,{Te  -»  r}  =  1  and  by  (2.12)  P(ot^(A)  — ►  1.  it 
suffices  to  show 

(2-14)  {1aW*(T)/(1  -  T);  r  >  0} 

is  uniformly  integrable. 

Let  t1  =  inf{t  :  W2(t)  >  |yi|/2}.  For  all  sufficiently  small  e  A  C  {T  <  r  /.  It  is  easy 
to  see  that  [Wi(t)  —  xo  —  (xx  —  x0)t]2/(l  —  <)2  —  t/(  1  — f),0  <  t  <  1,  is  a  martingale  and 
hence  [Wi(t)  —  i0  —  (xj  —  x<j)i]2/(l  —  t)2,0  <  t  <  1,  is  a  submartingale.  From  the  joint 
distribution  of  r'  and  Wj(r')  we  obtain 

£&!<;{1w'i(t')  _  10  _  (*> -  ^'iVd  -  02}  =  [»-7(i  -  r')i, 

which  is  finite  by  Lemma  1  below.  Also 

{["'■«  -  *0  -  (*i  -  *o)<l2/(i  -  «)2;  r'  >  (|  =  ((1  -  >0  —  0 

as  t  — ►  1,  again  by  Lemma  1.  It  follows  from  Doob’s  optional  sampling  theorem  that  on 
{T  <  t'} 

[W'.(T)  -  i0  -  (*i  -  io)T]2/(1  -  T)2 

<  E£yA  {[tV.d')  -  x„  -  (i,  -  x„)rf/(l  -  r')2|W'(f),l  <  t|. 

Hence  {l{r<r' )[Wi(T)  —  xo  —  (xx  —  xo)T]2/(l  —  T)2,e  >  0}  is  uniformly  integrable.  The 
uniform  integrability  of  (2.14)  follows  from  the  relation  A  C  {T  <  r'},  the  inequality 
(a  +  b )2  <  2(a2  +  b2),  and  Lemma  1. 
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Lemma  1  provides  justification  for  several  steps  in  the  preceding  argument  and  eval¬ 
uates  the  expectations  appearing  on  the  right  hand  side  of  (2.11).  It  will  be  convenient 
to  use  the  following  notation.  Let  W(t),0  <  <  <  oo,  denote  one  dimensional  Brownian 
motion  with  drift  \i  and  initial  value  W(0)  =  0.  We  write  PM  and  to  denote  dependence 
of  probabilities  and  expectations  on  fi.  For  b  >  0  let  r&  =  inf{<  :  W(t)  >  b},  where  it  is 
understood  that  inf  4>  =  -foo.  Let  =  Po(-|W(l)  =  f). 

Lemma  1.  For  0  <  t  <  1  the  P^(1)  density  function  of  Tb  is  given  by 

<2-15>  =  wrhw*  { 4 ^  - 6>  ( r^i)  ^ } 

For  £  >  b 

and 

E^[r2/(1  -  rb)\  =  6/tf  -  6)  -  »(- 0M0. 

Proof.  From  the  well  known  (and  easily  proved)  fact  that  the  P^  distribution  of 


W(-)  is  the  same  as  the  P$  distribution  of  (1  —  (-))W (  one  easdy  sees  that 

P^{r*  <t}  =  P^{(1  —  s)W[s/(l  —  s)]  >  b  for  some  s  <  t} 

=  Pt-b{n<t/(i-t)}. 

Equation  (2.15)  follows  by  differentiation  of  the  well  known  expression  for  the  last  proba¬ 
bility  (e.g.,  Siegmund,  1985,  (3.15)). 

From  (2.15)  one  obtains 

-n)]  =  »  jTV*/,(i  +»rV[({  - 

Writing  (1  +  s)~l  =  /0°°  interchanging  the  order  of  integration,  and  using  the 

well  known  equality 

e~atas~3^2(p(as~1^2  -  ^isx^2)ds  =  exp{-a[(2a  +  j*2)1/2  -  m]} 

i 

•> aa;v.%  a  .s.v  a  ••  .*•  '.  a  av.-  %.• ■  -  ^  v 


(e.g.,  Siegmund,  1985,  (3.16)  and  Problem  3.1),  one  obtains  the  given  expression  for 
E^\t* /(\  —  r*)].  A  similar  calculation  applies  to  ■E^(fi). 

Remarks,  (i)  As  observed  above,  the  boundary  of  7Z  can  be  defined  locally  near  (0. 0) 
by  a  function  y  =  /( x),  but  in  general  it  cannot  be  so  defined  globally.  However,  for  e 
sufficiently  small,  on  {  maxo<<<«  |W^(t)|  <  e1/4  }  T  defined  by  (2.1)  and  the  exit  time  from 
TZ  coincide,  so  there  is  no  loss  of  generality  restricting  attention  to  stopping  times  of  this 
form,  (ii)  In  higher  dimensions,  /"  in  (2.4)  becomes  the  Laplacian  A/,  and  x2(i  =  0,1) 
and  (xi  —  x0)2  become  Euclidean  distances  ||  Xi  ||2  (i  =  0, 1)  and  ||  x*  —  x0  ||2.  The  proof 
is  essentially  unchanged. 

In  Theorem  1  £o  and  are  at  a  distance  Ole1/2)  from  the  boundary  of  7Z  and  from 
each  other,  and  consequently  {T  <  e)  converges  to  a  limit  between  0  and  1.  Theorem 
2  is  concerned  with  the  case  that  £o  and  are  fixed  as  e  — ►  0,  so  {T  <  e)  — ►  0. 

As  above,  for  given  £<h£i  e  ^  suppose  there  exists  a  unique  point  on  7£,  the  sum  of 
whose  distances  from  £o  and  is  a  minimum,  and  consider  the  tangent-normal  coordinate 
system  through  this  point.  Let  £,  have  coordinates  (x,,y,)  (t  =  0, 1),  and  let  371  be  given 
by  the  graph  of  (x,  /(x))  in  some  neighborhood  of  (0,0),  so  /( 0)  =  /'( 0)  =  0. 


Theorem  2.  Assume  /  is  twice  continuously  differentiable,  y0yi  >  0,  and 
(2-16)  2y0yi/"(0)[l  +  (xi/yi)2]/|y0  +  y1|  >  -1. 


Let  T  =  inf{<  :  W(t)  e  371} .  Then  as  e  -*  0 


(2  17)  P(t)  {T  <  e\ _ _  — _ _ 

1  ■  ’  '  {1  +  2yoyi/"(0)(l  +  (-ri/yi  )3]/|yo  +  ytl}1^2 


exp(-2c  1 5/o£ri ) 


One  can  prove  Theorem  2  along  the  lines  of  the  proof  of  Theorem  1,  but  the  details 
are  rather  different.  To  keep  this  paper  to  a  reasonable  length  the  proof  has  been  omitted. 
An  example  comparing  the  numerical  accuracy  of  (2.17)  and  (2.4)  is  given  in  Section  4. 

An  interesting  case  which  fails  to  satisfy  the  conditions  of  Theorem  2  is  7?  a  disk  with 
£0  =  at  the  center.  In  this  case,  the  nearest  point  on  37Z  is  not  unique  and  (2.16)  is  not 
satisfied.  For  an  approximation  in  this  case,  which  leans  heavily  on  rotational  symmetry. 
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see  Siegmund  (1985,  Problem  11.1).  An  exact  expression  has  been  obtained  by  Kiefer 
(1959),  but  it  is  quite  complicated. 

A  related  but  somewhat  more  complicated  problem  than  that  discussed  in  Theorem 
1  is  to  approximate  the  joint  distribution  of  (X,  Wi(T)),  which  can  be  attacked  via  the 


characteristic  function 


(2.18) 


+  ,X2T/c];T  <  e). 


Expansion  of  (2.18)  to  the  precision  of  Theorem  1  seems  to  require  more  complicated 
calculations,  which  turn  out  to  be  very  similar  to  those  given  in  the  following  section  in 
order  to  obtain  the  term  of  order  e  in  the  expansion  of  {X  <  e}. 

It  seems  possible  to  obtain  the  results  of  this  section  by  the  methods  of  Jennen  and 
Lerche  (1981,1982),  but  the  computations  appear  to  be  somewhat  more  complicated.  If 
one  is  interested  in  the  joint  behavior  of  T  and  Wj(X),  their  method  might  turn  out  to  be 
the  simpler  one. 


3.  The  term  of  Order  e  and  Applications. 

Calculation  of  higher  order  terms  in  the  expansion  (2.4)  rapidly  becomes  very  com¬ 
plicated  in  detail.  In  this  section  we  see  what  is  involved  by  examining  the  term  of  order 
e.  (See  equation  (3.10).)  To  simplify  the  algebra  we  suppose  that  =  £i-  This  special 
case  sufficies  for  applications  to  the  problems  of  Kac  (1966)  and  Handelsman  and  Keller 
(1967),  which  are  discussed  below  (cf.  (3.15)).  We  proceed  informally  as  in  the  first  part 
of  the  proof  of  Theorem  1.  The  localization  and  uniform  integrability  arguments  necessary 
for  a  rigorous  proof  are  similar  to  those  in  Theorem  1  and  have  been  omitted. 

Let  £0  =  £i •  In  the  notation  of  Section  2  for  the  standardized  problem  on  the  time 
interval  [0,1],  (2.6)  becomes 

<31>  PmAT  <  l}^p(2Vo)  =  |exp  ^ - JT/TfYZTf) - j  /  ' 

where  £o  =  (0,yo),£o  =  (0,  |yo|),  and  T  =  Te  is  defined  by  (2.3).  Assuming  that  /  is  three 
times  continuously  differentiable,  we  have 

E~1/2f(e,/2x)  =  e'/2f"(  0)x2/2  +  ef'"(  0)x3/6  +  o(s); 


and  hence  the  right  hand  side  of  (3.1)  becomes 


(3.2)  i  -  -  t)] 

-  i£|yol/"'(0)£«;e.(^(T)/(l  - T)]  +  [W?(T)/(1  - T)2] 

+  o(£). 

Until  further  notice,  we  shall  write  P  and  E  for  PiX\,  and  Eix\, .  Recall  the  definition 

vO»So  sO»Qo 

of  r  given  in  (2.9)  and  note  that  the  conditional  distribution  of  W\(t)  given  r  is  normal 
with  mean  0  and  variance  r(l  —  r).  Since  P{Te  — ►  r}  =  1  (e  — ♦  0),  we  have 

(3.3)  E[W*(T)/(  1  -  T)]  -  0 
and 

(3.4)  E[W}(T)/(  1  -  T)2]  -►  3 Et2. 

Also 

(3.5)  £[tV?(T)/(l  -  T)]  =  Er  +  {E[WftT)l(\  -  T))  -  E(H?( r)/(l  -  r)]}; 

and  the  final  contribution  to  the  term  of  order  e  in  (3.2)  comes  from  the  difference  on  the 
right  hand  side  of  (3.5),  which  is  itself  of  order  e1/2. 

First  suppose  that  /"(0)  <  0  and  to  simplify  some  details  that  /(x)  <  0  for  all  x.  The 
case  /"( 0)  >  0  involves  a  similar  argument  with  slightly  more  complicated  calculations. 
Let  Tt  denote  the  <r-algebra  generated  by  W(.s),  s  <  t.  Since  T  <  r,  we  have 

(3.6)  £«( r)/(l  -  T )]  =  £[(IV,(r)  +  W\(t)  -  tV,(T))2/(l  -  r)l 

=  E  {  22^1  +  ^^E[W,( r)  -  Wi{T)\Tt,  r]  +  (1  -  r)-‘  E[(W,(t)  -  WMT))2  \FT.  r\ J 

Conditional  on  Tt  and  r,  Wi(r)  —  W\(T)  is  normally  distributed  with  mean  —  Wx(T)(t  - 
T)/(  1  —  T)  and  variance  (r  —  T)(l  —  r)/(  1  —  T).  Hence  after  some  algebra  one  obtains 

(3.7)  E[W*(T)/(  1  -  T)\  -  E[W*(t)I(\  -  r)] 

=  £[(1  -  T)~2W2(T)E(t  -  T\Tt))  -  E{(  1  -  T)~l E(r  -  T\Tt)}. 
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Doob’s  optional  sampling  theorem  yields 


£■{(1  -  t)-'[W2(t)  -  W,(T)  -  (r  -  T)(|y„|  -  W,(T))/(  1  -  T)\ \?T)  =  0 


(cf.  Siegmund,  1985,  Problem  3.12),  and  hence  with  probability  one  as  e 

^■/al/"(Q)|tv?(r)(i-r) 


(3.8) 


E(r  -  T\?t) 


|y„l  +  it‘/J|/"(o)|w?(r) 


-=1/2|yor1|/"(0)|W,2(r)(l  -  r). 


1-1 1  e» 


Substitution  of  (3.8)  into  (3.7)  yields 

(3.9)  E{W}(T)/(1  -  T)|  -  £[W\2( r)/(l  -  r)] 

~5*'/,|!»r,ir(0)l{£«(r)/(l-r)|-£:W'1V)} 

=  5«l/J|y.r,IA0)H3E[r*(l  -  r))  -  £(r(l  -  r)|}. 

From  (3.1)-(3.5)  and  (3.9)  we  finally  obtain 

(3.10)  P™,{T  <  1}  =  «p(-2y’){l  -  e,/,|y.ir(«)^!e;(r) 

+  iE(/"(0)|I£<<;,<,(3rI(l  -  r)  -  r(l  -  r)  +  3y0V]  +  o(e)|, 

where  T  is  defined  by  (2.3),  r  by  (2.9)  £o  =  (0,  yo),  and  £o  =  (0, 1 y o  I ) - 

The  expansion  (3.10)  also  holds  when  /"( 0)  >  0.  In  this  case  r  <  T,  so  (3.6)  must  be 
replaced  by 


E\W*(T)/(1  -  T)]  =  E\Wl(r)/(\  -  T)| 

W 


+2£  ( Wl(r)E  [^(X)_-A(r)^]  |+E|£  } 


By  optional  sampling 

'Wx{T)-Wx{t) 

1  -  T 


I  Ft 


=  _(1  _  t)-'Wx(t)E[(T  -  T  )/( 1  -  T) \Tr), 


and  it  may  be  shown  that 


E{[W,(T)  -  W,(r)]V(  1  -  T) \Fr)  ~  £[(T  -  r)/(l  -  T)|^r], 

Hence  in  place  of  the  equality  (3.7)  one  obtains 

EW(t)H  1  -  r)]  -  E[Wl(T)H  1  -  T)1 
~  £|(1  -  t)~2W%(t)E(T  -  r\?r)\  -  £[fl  -  t)-'E(T  -  r\Tr)\. 

A  result  similar  to  (3.8)  holds  for  E(T  —  rj^>),  and  the  rest  follows  as  before. 

Remarks,  (i)  By  the  method  of  Lemma  1  one  can  evaluate  the  moments  appearing 
on  the  right  hand  side  of  (3.1).  However,  for  the  applications  given  below,  which  in  effect 
involve  an  integration  of  (3.10)  over  the  computations  are  considerably  simpler  if  one 
interchanges  the  order  of  the  two  integrations  and  integrates  over  £0  first,  (ii)  In  higher 
dimensions  the  relation  of  &R.  to  its  tangent  planes  can  be  more  complicated  than  in  two 
dimensions.  In  general,  one  must  condition  on  Tt^t  and  consider  the  two  cases  {T  <  r} 
and  {T  >  r}.  Whereas  the  term  of  order  e1/2  involves  only  the  Laplacian  of  /,  i.e.,  the 
mean  curvature  of  &R,  the  term  of  order  e  involves  mixed  partials  as  well.  For  the  problem 
studied  by  Handelsman  and  Keller  (1966),  where  H  is  the  region  exterior  to  a  sphere  in 
2R3  one  does  not  encounter  these  complications. 

Now  let  T  denote  the  first  exit  time  of  W  from  R,  and  for  £0,  £i  e  H  define  p(t,  £0,6) 
by 

p(Uo,6)#i  =  Pr (T  >  t,W(t)  e  d^lWiO)  =  to)- 

Observe  that 

(3.U)  =  (2-rt)-"2[l  -  PW{0{T  <  <}]. 

In  order  to  study  (1.2)  in  a  bounded  region  7 Z  in  2R 2  Kac  (1966)  uses  the  representation 

^exp(-Afct)  =  JJ*p(t,(o,to)dt  o, 

which  by  (3.11)  equals 

(27rt)-,(|7Z|  -  f <  t)dU 


(3.12) 


where  \7l\  denots  the  area  of  71.  Handelsman  and  Keller  (1966)  are  interested  in  1R3  and 
the  integral 

JJJji1  ~  (2xt)3/2p(<,fo,6>Mo, 

which  by  (3.11)  equals 

(313)  ////“(T<‘R"' 

In  order  to  analyze  the  integral  in  (3.12)  it  is  convenient  to  make  a  change  of  variables 
(cf.,  Pleijel,  1954)  to  obtain 

(3.14)  f f  p£(q{T  <  t}d( o  =  /  /  P5(.°{.{r  <  1}[1  -  \yo\c{a)]d\yo\da  +  0(e-2<>), 

J  J  H  J  dTt,  J  0 

where  a  denotes  arc  length  on  dH,c(‘)  is  the  curvature  of  dll,  and  fo  has  coordinates 
(0,  i/o )  in  the  tangent-normal  coordinate  system  with  its  origin  at  the  point  a  of  dll.  so 
|yo|  is  the  distance  from  £0  to  dR  and  c(<x)  =  —  /M(0). 

Keeping  (2.2)  and  (2.15)  in  mind,  one  can  substitute  (3.10)  into  (3.14),  integrate  with 
respect  to  |yo|,  then  with  respect  to  <r,  and  refer  to  the  Gauss-Bonnet  theorem  as  indicated 
by  Kac  to  obtain 

(3.15)  £  exp(-A kt)  =  (27Tt)"1  \H\  -  [4(2ir<)1/2]"1  \dfc\ 

+(1  -  h)/6  +  2“8(27t)-1/2  t1'2  +  o(t1'2), 

where  \&R.\  is  the  length  of  dR  and  h  is  the  number  of  holes  in  71. 

Since  (3.15)  involves  integration  of  (3.10),  some  additional  justification  is  required  to 
claim  that  (3.15)  has  been  proved  rigorously.  This  seems  a  straightforward,  albeit  rather 
technical  matter.  Since  it  does  not  appear  to  add  significant  insight,  the  details  are  omitted. 

The  expansion  (3.15)  agrees  with  those  given  by  Stewartson  and  Waechter  (1971)  and 
Smith  (1981),  both  of  whom  used  analytical  methods  and  obtained  additional  terms.  The 
term  of  order  t1!2  disagrees  with  that  given  by  Louchaxd  (1968),  whose  argument  appears 
to  contain  an  improper  use  of  the  Markov  property. 

A  similar  computation  yields  the  expansion  of  Handelsman  and  Keller  (1966). 


4.  Monte  Carlo  Methods. 

Again  let  W(t )  denote  Brownian  motion  starting  from  some  point  inside  a  region  1Z  in 
d-dimensional  Euclidean  space,  and  let  T  =  inf{t  :  W(t)  4  71).  In  this  section  we  consider 
the  problem  of  estimating  by  Monte  Carlo  methods  probabilities  like 

(4.1)  Pr{T  <  <}. 

The  same  ideas  axe  applicable  to  substantially  more  complicated  first  passage  distributions. 

An  obvious  procedure  to  estimate  (4.1)  is  to  partition  the  time  interval  [0,  t]  by  the 
points  ti  =  iff,  i  =  0, 1, . . . ,  m,  where  e  =  </m,  generate  N  realizations  of  the  discrete  time 
random  walk  W(fj),  *  =  0, 1, ...  ,m,  and  estimate  (4.1)  by  the  relative  frequency  among 
the  N  realizations  that  W(<i)  4  ft  for  some  1  <  t  <  m. 

The  standard  deviation  of  this  estimator  is  of  order  TV-1/2.  Its  bias  equals  the  dif¬ 
ference  between  Pr{W(tj)  4  ft  for  some  1  <  i  <  m}  and  (4.1),  which  presumably  is  of 
order  e1/2  (cf.  Nagaev,  1970,  Siegmund,  1985,  Chapter  X,  Hogan,  1984).  Thus  the  bias 
is  of  the  same  order  as  the  sampling  error  unless  e  is  small  compared  to  IV-1 .  Since  N 
may  be  in  the  thousands,  it  is  often  computationally  unfeasible  to  achieve  a  satisfactory 
estimate  by  the  obvious  device  of  making  e  extremely  small. 

The  procedure  we  propose  to  study  is  the  following.  Having  generated  the  partial  real¬ 
ization  W(tQ), ... ,  W(ti)  and  decided  that  T  >  tj,  generate  W(ti+i  )•  If  W  (<«-t-i )  4  ft  decide 
T  <  <j+i  <  t.  If  W(U+i)  e  ft  decide  T  <  tj+i  <  t  with  probability  p[W(*i), IV(<j+i),ff], 
where  p(^0i^i>c)  is  a  suitable  approximation  for  {T  <  ff}  obtained  from  Theorem  1. 

As  a  first  example  we  consider  a  rather  complicated,  but  linear  problem.  In  this  case 
there  is  no  question  how  carefully  one  should  approximate  P\[\x  {T  <  e},  which  can  be 
evaluated  exactly;  and  we  see  that  a  striking  improvement  in  the  accuracy  of  the  naive 
estimator  is  possible. 

The  example  concerns 

(4.2)  -  W(,»  >  6), 

where  W  is  one  dimensional  Brownian  motion  and  b  is  substantially  larger  than  max(0.  £). 
so  the  probability  (4.2)  is  small.  This  probability  has  arisen  in  the  unrelated  problems  of 
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Levin  and  Kline  (1985)  and  Adler  and  Brown  (1986).  Hogan  and  Siegmund  (1986)  show 
that  if  b  — ►  oo  and  £  varies  with  b  in  such  a  way  that  £/6  is  a  fixed  real  number  less  than 
1,  then  (4.2)  equals 

(4.3)  [2(26  -  {)(6  -  {)  +  1  +  o(l)]  exp[-26(6  -  {)]. 

In  order  to  check  the  accuracy  of  the  approximation  (4.3)  by  simulation,  we  first 
generate  a  Pq  ^  realization  of  the  discrete  skeleton  W(i/m),  i  =  0, 1, . . . ,  m  and  locate  the 
points  0  <  ui  <  1/2  <  m  which  satisfy 

W(v2/m)  —  W{v\/m)  =  max  [W(j /m)  —  W(i/m)\. 

0  <i<j<m 

Then  for  each  [(t/i  4-  1/2)/ 2]  <  i  <  m  we  generate  the  maximum  of  W(<)  for  i/m  <  t  < 
(i  +  1  )/m,  which  conditional  on  W(i/m)  —  £0,  ^{(i  +  1)/”*}  =  £1  has  the  distribution 

P\Ztx  =  exp[-2m(x  -  £o)(x  -  6)] 

forx  >  max(^0,^i)-  Similarly  we  generate  minima  over  [i/m,  (i+ 1  )/m],  i  =  0,1 . [(^i  + 

u2)/2]  —  1.  Putting  /*  =  m_1[(/ix  +  ^2)/2],  we  use 

(4.4)  max  W(t)  —  min  W(<) 

r<i<i  o<i<c* 

as  a  surrogate  for  the  desired 

max  [W(<)  —  W(s)l, 

°<*<‘ 

for  which  it  is  in  fact  a  lower  bound.  Presumably  the  discrepancy  between  these  two 

quantities  is  unimportant  in  the  cases  of  primary  interest,  when  (4.2)  is  small.  As  a  check 

one  might  compare  results  for  two  different  values  of  m,  or  alternatively  perform  a  second 

experiment  with  the  maxima  and  minima  taken  over  overlapping  sets  of  intervals,  say 

V\  <i  <m  and  0  <  t  <  i/2,  which  in  all  but  very  few  cases  would  yield  an  upper  bound. 

Table  1  gives  the  results  of  an  experiment  with  N  =  9999  repetitions.  The  first  Monte 

Carlo  estimate  reported  in  each  row  is  the  relative  frequency  of  the  event  {maxo<l<;<m[H  t  J , 
« 

W(i/m)\  >  b };  the  second  is  the  modified  estimate  described  above.  The  final  entry  in 
each  row  gives  the  approximation  (4.3).  There  is  wide  disparity  between  the  discrete 
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skeleton  estimator  and  the  modified  estimator,  and  except  for  two  cases  there  is  excellent 
agreement  between  the  latter  estimator  and  the  theoretical  approximation.  Those  two  dis¬ 
crepancies  both  involve  large  probabilities,  where  (4.3)  is  not  expected  to  provide  a  good 
approximation. 

For  a  second  example  we  consider  the  first  hitting  time  of  a  sphere  of  radius  r  by  a 
three  dimensional  Brownian  motion  starting  outside  the  sphere.  This  problem  is  surrogate 
for  a  much  more  elaborate  problem  in  physical  chemistry  (N.  J.B.  Green,  personal  commu¬ 
nication).  In  that  problem  n  independent  spheres  of  radius  r'  follow  independent  Brownian 
paths,  annihilating  each  other  if  they  collide.  Our  problem  is  the  spescial  case  n  =  2  and 
r  =  2 r'.  Moreover,  if  one  simulates  a  sequence  of  snapshots  of  the  configurations  of  the  n 
spheres  at  times  ie  ( *  =  0, 1, . . . , .),  it  seems  plausible  that  one  can  bridge  the  short  gap 
from  ie  to  (i  +  l)c  by  considering  each  pair  of  spheres  in  isolation  from  the  others,  and 
hence  the  case  n  =  2  may  be  useful  preparation  for  other  cases. 

In  order  to  implement  the  proposed  algorithm  one  must  choose  an  approximation  for 

(“•s>  ' '(1  (  mm  ||  tV(f)ll<rl 

Since  the  approximation  of  Theorem  1  requires  some  numerical  computation  to  determine 
the  point  on  the  surface  of  the  sphere  the  sum  of  whose  distances  from  £o  and  £i  is  a 
minimum,  it  seems  reasonable  to  try  first  the  simpler  approximation  which  treats  the 
surface  of  the  sphere  as  a  plane  and  approximates  (4.5)  by  the  very  simple 

(4-6)  exp[-2(||  6,  11  ~r)( ||  6  ||  -r)/e]. 

Under  the  asymptotic  scaling  of  Theorem  1,  (4.6)  and  (2.4)  have  the  same  limits,  but  differ 
at  the  term  of  order  e1^7. 

As  indicated  above,  there  is  reason  to  think  that  the  bias  of  a  direct  frequency  count 
from  the  discrete  skeleton  is  of  order  e1^7  and  hence  requires  one  to  take  e  =*  JV-1.  the 
number  of  repetitions  of  the  experiment.  Although  we  have  no  idea  how  to  prove  it,  we 
believe  that  use  of  (4.6)  to  bridge  the  gaps  in  the  discrete  skeleton  results  in  an  estimator 
whose  bias  is  of  order  e,  and  use  of  the  better  approximation  of  Theorem  1  reduces  this 
bias  to  order  e3/2  .  If  these  conjectures  are  true,  then  even  if  the  simple  approximation 


(4.6)  is  used,  a  reasonable  magnitude  for  e  is  about  iV-1/2,  which  represents  a  considerable 
improvement. 

Table  2  reports  the  results  of  a  simulation  experiment  to  estimate  (4.1)  for  the  hitting 
time  T  of  a  sphere  of  radius  r  by  a  three  dimensional  Brownian  motion  starting  from 
the  point  (x,  0, 0).  The  interval  width  of  the  discrete  skeleton  is  £.  The  Monte  Carlo 
estimators  are  based  on  N  =  9999  repetitions  of  the  method  described  above  with  the 
simple  approximation  (4.6).  The  exact  probabilities, 

Pr(T  <  t)  =  2x-1r{l  -  *[(x  -  r)/t1/2]}, 


are  also  included.  The  approximations  are  quite  accurate,  in  fact  more  accurate  than  one 
would  expect  from  the  conjectured  size  of  the  bias. 

Our  final  example  is  designed  to  discover  whether  we  can  gain  anything  by  using  the 
presumably  better  approximation  provded  by  Theorem  1  to  bridge  the  gap  between  H'(  t, ) 
and  W(/j+i).  To  simplify  the  programming  problem  of  determining  the  local  tangent - 
normal  coordinate  systems  associated  with  the  points  W(ti)  =  and  W(tI+1)  =  .  we 

take  d  —  2  and  T  the  first  time  W(t)  is  within  a  circle  of  radius  r  centered  at  the  origin. 
In  order  to  investigate  at  the  same  time  the  numerical  accuracy  of  the  approximations  of 
Theorems  1  and  2,  we  study  the  conditional  probabililty,  Pfl\0{T  <  1}. 

For  V.  the  region  exterior  to  a  circle  of  radius  r, /"( 0)  =  r-1.  It  seems  plausible  that 
one  might  improve  the  numerical  accuracy  of  (2.4)  slightly  by  including  the  term  of  order 
in  the  exponent  to  obtain 


<4-7) 


|  -  2t_1  lyoVi |  -  r  1  t  1/J|sro»il 


*[g~1/2(yo  +  vi )] 
v>k-1/2(yo  +  yi)] 


[(1  -  r  ‘(xo  -  *i)2]  +e"1(io|yi|  +*?|yol) 


Table  3  contains  the  results  of  an  experiment  with  r  =  1  and  N  =  9999  repetitions. 
There  is  no  evidence  that  the  more  complicated  approximation  yields  more  accurate  Monte 
Carlo  estimates.  Since  the  additional  numerical  computation  required  to  evaluate  (4.7) 
after  each  new  observation  more  than  doubles  the  total  computational  effort,  there  does 


not  seem  to  be  any  reason  to  use  this  approximation,  at  least  for  the  relatively  simple 
problem  considered  here. 

As  expected,  the  analytical  approximation  provided  by  Theorem  1  is  better  than  that 
of  Theorem  2  when  P^(0  {T  <  1}  is  large  or  moderate,  and  the  converse  is  true  when  this 
probability  is  small.  Overall,  the  better  of  the  two  approximations  is  reasonably  good. 
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Table  1 


Evaluation  of  {maxo<4«<i  [W(t)  —  W(s))  >  b} 


Table  2 

First  Passage  to  a  Sphere 


r 

X 

£ 

t 

Pr(T  <  t) 

f 

Pr(T  <  t') 

Monte  Carlo 

Exact 

Monte  Carlo 

Exact 

2 

3 

.50 

1 

.2110 

.2115 

9 

.5016 

.4926 

2 

3 

.25 

1 

.2162 

.2115 

9 

.4996 

.4926 

4 

6 

.50 

2 

.0986 

.1049 

25 

.4495 

.4594 

4 

7 

.25 

2 

.1004 

.1049 

25 

.4597 

.4594 
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Table  3 

Approximations  for  P^*^  {mino<*<i  ||  W(t)  ||<  1} 

||  £o  ||  Monte  Carlo  Estimates  Analytic  Approximations 

Using  (4.6)  Using  (4.7)  (4.7)  (2.17) 


£  =  .1  £  =  .05  £  =  .1  £  =  .05 


1.2 

.866 

.891 

.856 

.889 

.889 

.343 

1.5 

.536 

.535 

.524 

.538 

.515 

.495 

2.0 

.097 

.106 

.095 

.102 

.089 

.090 
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